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Abstract 

The  dynamics  of  solid  oxide  fuel  cell  (SOFC)  operation  have  been  considered  previously,  but  mainly  through  the  use  of  one-dimensional 
codes  applied  to  co-flow  fuel  cell  systems.  In  this  paper  several  geometries  are  considered,  including  cross-flow,  co-flow,  and  counter-flow. 
The  details  of  the  model  are  provided,  and  the  model  is  compared  with  some  initial  experimental  data.  For  parameters  typical  of  SOFC 
operation,  a  variety  of  transient  cases  are  investigated,  including  representative  load  increase  and  decrease  and  system  shutdown.  Of  particular 
note  for  large  load  decrease  conditions  (e.g.,  shutdown)  is  the  occurrence  of  reverse  current  over  significant  portions  of  the  cell,  starting  from 
the  moment  of  load  loss  up  to  the  point  where  equilibrated  conditions  again  provide  positive  current.  Consideration  is  given  as  to  when  such 
reverse  current  conditions  might  most  significantly  impact  the  reliability  of  the  cell. 
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1.  Introduction 

Development  of  fuel  cell  systems  is  being  heavily  pursued 
to  help  provide  society’s  future  energy  needs.  Thus  far  there 
have  been  many  demonstrations  for  a  variety  of  fuel  cell  tech¬ 
nologies,  and  some  units  are  now  commercial.  Currently  the 
U.S.  Department  of  Energy  is  sponsoring  the  development  of 
solid  oxide  fuel  cell  technology  through  its  Solid  State  Energy 
Conversion  Alliance  (SEC  A)  Program.  This  program  focuses 
on  solid  oxide  fuel  cell  (SOFC)  technology  with  power  output 
between  3  and  10  kW  and  having  designs  suitable  for  low  cost 
mass  production  [1].  Enabling  mass  production  capability  for 
fuel  cell  technology  will  help  to  reduce  their  present  high 
capital  cost  (~$  5000  kW-1)  to  less  than  $  400  kW-1,  mak¬ 
ing  them  economically  competitive  with  other  distributed 
power  generation  methods.  In  addition,  because  of  their 
superior  emissions  performance,  the  environmental  permit 
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process  is  expected  to  be  relatively  easy.  As  a  result,  they 
can  be  readily  applied  for  such  applications  as  grid- support 
[2]. 

Most  of  the  focus  in  the  development  of  fuel  cell  technol¬ 
ogy  has  been  on  achieving  good  steady  state  performance. 
Good  performance  includes  not  only  good  fuel  to  electric 
conversion  efficiency,  but  also  long  lifetime  (e.g.,  >40,000  h 
for  stationary  applications).  While  good  steady  state  per¬ 
formance  is  critical  for  the  commercialization  of  fuel  cell 
technology,  so  also  is  dynamic  performance.  Given  that  real- 
world  applications  will  incur  frequent  load  changes,  partic¬ 
ularly  at  the  low  power  application  range  (^3-30  kW),  sys¬ 
tems  targeting  such  applications  will  likely  need  to  have  load¬ 
following  capability  if  they  are  to  be  commercially  viable. 
Numerous  investigators  have  examined  the  overall  perfor¬ 
mance  of  fuel  cell  systems,  and  some  have  examined  sys¬ 
tem  dynamics  due  to  load  changes,  e.g.  [3-7].  However,  at 
least  for  SOFC  technology,  very  little  has  yet  been  done  to 
examine  what  occurs  within  the  cell/stack  due  to  load  vari¬ 
ations  [8-10].  One  reason  for  the  lack  of  work  in  this  area 
relates  to  the  difficulty  of  experimentally  accessing  the  in- 
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Nomenclature 

A  area  at  the  inlet/outlet  of  the  control  volume 

(m2) 

Aa  electrochemical  active  area  of  the  control  vol¬ 
ume  (m2) 

Ci  molar  concentration  of  species  i  (kg  mol  m-3) 
Qo  molar  concentration  of  species  i  in  the 

freestream 

E  specific  energy  (J  kg-1);  Nernst  voltage  (V) 

E  rate  of  energy  across  control  surface  due  to 
electrochemical  reaction 
Lceii  operating  cell  voltage  (V) 

F  Faraday’s  constant 

(96,485  C  g  mol-1  electrons) 

Fs  control  volume  surface  frictional  force  (N) 

G  Gibbs  free  energy 

h  flow  channel  height 

i  current  density  (A  cm-2) 

i0  exchange  current  density  (A  cm-2) 

iL  diffusion  limiting  current  density 

La  voltage  loss  due  to  electrochemical  reaction 

(V) 

Lc  voltage  loss  due  to  concentration  losses  result¬ 
ing  from  diffusion  of  reactants  (V) 

Lr  voltage  loss  due  to  ionic/electronic  resistive 
media  (V) 

M  momentum  within  control  volume  (kg  m2  s- 1 ) 
M  rate  of  momentum  flow  across  a  control  surface 

(kg  ms-2) 

n  number  of  participating  electrons  in  fundamen¬ 

tal  electrochemical  reaction 
Ni  molar  flow  rate  of  species  i  (in/out  of  control 
volume)  (kg  mol  s  - 1 ) 

P  pressure  at  the  inlet/outlet  of  the  control  vol¬ 

ume  (Pa) 

P'  net  rate  of  flow  energy  added  at  control  volume 
inlet  and  outlet.  Zero  for  exit  control  volume 
where  enthalpy  is  directly  used  in  the  conser¬ 
vation  equation 

gs  net  boundary  heat  addition  (W) 

R  net  resistance  of  cell  ( =R0  +RC)  (Q  m2) 

R0  ohmic  resistance  of  the  cell  (Q  m2) 

Rc  contact  resistance  between  cell  conductive  lay¬ 

ers  (Q  m2) 

Ri  molar  source  of  species  i  due  to  electrochemi¬ 

cal  reaction  (kg  mol  s-1) 

Ru  universal  gas  constant 

S  supply  flow  rate 

AS  entropy  generated  per  mole  of  reactant 

(Jgmol-1  K-1) 
t  material  thickness 

T  temperature  (K) 

v  velocity  (ms-1) 

V_ _ finite  control  volume  (m3) _ 


~w  mass  flow  (kg  s  1 ) 

y  mole  fraction 

Greek  symbols 

a  electrochemical  transfer  coefficient  [1] 

0  electrochemical  overpotential  (V) 

p  density  (kg  m- 3 ) 

Subscripts  and  superscripts 
a  fuel  electrode  (anode  for  i  >  0  A  cm-2) 

Air  air 

c  air  electrode  (cathode  for  i  >  0  A  cm-2) 

e  electrode 

H2  hydrogen 

i  species  index 


ternal  features  of  a  cell/stack  without  modifying  the  flow 
and  temperature  which  are  critical  to  determining  the  cell 
performance. 

At  the  National  Energy  Technology  Laboratory  (NETL), 
work  is  underway  to  examine  cell  and  stack  dynamic  perfor¬ 
mance.  Both  models  and  experimental  work  are  being  pur¬ 
sued.  Experimentally,  we  have  been  studying  the  ‘random’ 
fluctuations  that  occur  in  PEM  fuel  cells  when  heavily  loaded 
and  methods  for  reducing  such  fluctuations  in  order  to  extend 
the  operating  range  of  a  fuel  cell  [11].  This  paper  reports 
on  modeling  efforts  related  to  SOFC  cell/stack  performance 
under  transient  loads.  Three  planar  geometries  are  consid¬ 
ered:  co-flow,  cross-flow,  and  counter-flow  (CoF,  CrF,  CtF, 
respectively).  In  the  following  section,  a  detailed  description 
is  given  of  the  model  tool  developed  at  NETL  for  dynamic 
fuel  cell  analysis.  One  previous  application  of  this  tool  has 
been  in  the  study  of  the  dynamic  performance  of  hybrid  sys¬ 
tems  [3,12].  In  Section  4,  discussion  of  how  the  model  is 
applied  to  investigate  the  load  transient  characteristics  of  an 
SOFC  is  made.  Next,  results  are  given,  and  following  which 
a  detailed  analysis  of  the  results  is  presented.  Finally,  conclu¬ 
sions  from  the  present  work  are  made,  along  with  suggestions 
for  future  work. 


2.  Model  description 

In  the  present  discussion,  a  generic  bi-polar,  planar  fuel 
cell  is  examined,  see  Fig.  1.  Straight,  single  pass  flow  chan¬ 
nels  are  assumed  to  distribute  fuel  and  oxidant  to  the  fuel 
electrode  and  air  electrode  sides  of  the  cell,  respectively.  (We 
avoid  the  terms  anode  and  cathode  for  reasons  that  will  be 
understood  later.)  As  described  in  detail  below,  various  sub¬ 
models  can  be  employed  for  analyzing  the  conditions  within 
these  channels  depending  on  the  requirements  of  the  model 
application  to  resolve  various  time  scales.  Regardless  of  flow 
geometry,  CoF,  CrF,  or  CtF,  these  flow  channels  can  be  an¬ 
alyzed  using  one-dimensional  flow  calculations  to  provide 
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Fig.  1.  Planar  cross-flow  cell  geometry. 


useful  design  information.  While  CoF  and  CtF  geometries 
can  be  analyzed  using  single  channels  for  the  fuel  and  air, 
the  cross-flow  model  must  be  analyzed  using  an  assembly 
of  one-dimensional  flow  channels  [13].  The  justification  for 
using  one-dimensional  flow  calculations  can  also  be  found 
from  an  analysis  of  the  characteristic  time  scales  present  in  a 
fuel  cell. 


2.7.  Characteristic  time  scales 

A  first-order  analysis  provides  estimates  for  the  various 
time  scales  occurring  in  SOFC  systems.  Results  from  such 
an  analysis  are  given  in  Table  1 .  The  analysis  assumes  an  an¬ 
ode  supported  SOFC  fuel  cell  having  cell  parameters  shown 
in  the  last  column  of  Table  3.  The  first  column  in  Table  1 
is  the  name  of  a  particular  process  time  scale.  Associated 
with  each  process  is  its  rate  and  capacity  parameters,  and 
equations  are  provided  in  the  table  for  each  value.  Here, 
p  is  the  density  (mass  or  mol),  Ax  the  flow  passage  cross- 
section,  V  the  flow  passage  volume,  Vc  the  cell  voltage,  Ac 
the  active  area,  D  the  molecular  diffusion  coefficient,  Cp 
the  specific  heat,  C  the  electric  or  thermal  capacity,  t  the 
electrode  thickness,  L  the  cell  length,  and  V  the  gas  veloc¬ 
ity.  The  ratio  of  the  capacity  and  rate  values  provide  the 
first-order  estimate  of  the  time  scale  for  each  process.  Be¬ 
cause  some  transport  processes  depend  on  a  characteristic 
length  scale  (e.g.,  diffusion),  and  because  there  are  several 
distinct  characteristic  length  scales  within  a  fuel  cell  (elec¬ 
trode  thickness,  cell  length,  etc.),  numerous  transient  time 
scales  are  present  for  the  same  fundamental  transport  mech¬ 
anism. 


Table  1 


Characteristic  transient  time  scales  of  a  planar  SOFC 


Time  scale  name 

Process  rate  parameter 

Capacity  parameter 

Time 

Name 

Value 

Name 

Value 

scale,  t  (s) 

A.  Cell  charging  time 

Current  density  (I/Ac) 

104  Am-2 

Double-layer  capacity 
(CV c) 

10-1  Cm-2 

nr5 

B.  Cathode  gas  electrode 
mass  diffusion  time 

Molecular  diffusion 
rate  ( pDAJt ) 

10_1  mol  s-1 

Cathode  electrode  mole 
capacity  (ptAc) 

1CT6  mol 

nr5 

C.  Anode  gas  electrode 
mass  diffusion  time 

Molecular  diffusion 
rate  ( pDAJt ) 

10-1  mol  s-1 

Anode  electrode  mole 
capacity  (ptAc) 

1CT4  mol 

1(T3 

D.  Cathode  electrode 

thermal  diffusion  time 

Thermal  diffusion  rate 
(kAJf) 

103  Js"1  K-1 

Cathode  electrode 
thermal  capacity  (pCptAc) 

10°  JK"1 

1(T3 

E.  Convective  transport 
time 

Specific  gas  volume 
flow  rate  ( V/Ax ) 

10°  ms-1 

Specific  gas  volume 
(V/Ax) 

10“ 1  m 

0"1 

F.  Cell  Reactant 
consumption  time 

Faradaic  mole  density 
(I/(nFAc)) 

10~2  mol  (s  m2)-1 

Specific  mole  density 
(ph/MW) 

1CT2  molm-2 

10° 

G.Anode  gas  streamwise 
mass  diffusion  time 

Molecular  diffusion 
rate  ( pDAJL ) 

KT5  mol  s_1 

Anode  cell  mole  capacity 
(pLAx) 

KT4  mol 

101 

H.  Anode  gas  streamwise 
thermal  diffusion  time 

Thermal  diffusion  rate 
( kAx/L ) 

10-4Js-1  K-1 

Anode- gas  thermal 
capacity  (pCpLAx) 

10-3  JK-1 

101 

I.  Cathode  gas  streamwise 
thermal  diffusion  time 

Thermal  diffusion  rate 
( kAx/L ) 

10-4  J  s-1  K"1 

Cathode  gas  thermal 
capacity  (pCpLAx) 

10-3  JK-1 

101 

J.  Anode  electrode 

thermal  diffusion  time 

Thermal  diffusion  rate 
(kAJf) 

101  J  s-1  K-1 

Anode  electrode  thermal 
capacity  (pCptAc) 

101  JK"1 

101 

K.  Cathode  gas 
streamwise  mass 

Molecular  diffusion 
rate  ( pDAJL ) 

10~5  mol  s-1 

Cathode  cell  mole 
capacity  (pLAx) 

1CT4  mol 

101 

diffusion  time 

L.  Cell  heating  time 

Gas  heat  transfer  rate 

(Me) 

10-2  J  s-1  K-1 

Cell  heat  capacity  (pCtAJ) 

101  JK-1 

103 

M.  Anode  electrode 

streamwise  thermal 

Thermal  diffusion  rate 
(kAJL) 

10-2  J  s-1  K"1 

Anode  electrode  thermal 
capacity  (pCLAx) 

101  JK"1 

103 

diffusion  time 

N.  Cathode  electrode 

streamwise  thermal 

Thermal  diffusion  rate 
(kAJL) 

10-4  J  s-1  K“' 

Cathode  electrode  thermal 
capacity  (pCLAx) 

10°  JK-1 

104 

diffusion  time 
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The  results  of  Table  1  show  a  wide  range  of  time  scales 
(10-5  to  104  s)  spanning  over  nine  orders  of  magnitude.  These 
results  can  be  helpful  to  guide  engineering  analysis  by  show¬ 
ing  how  transport  models  can  be  simplified,  and  yet  still  ac¬ 
curately  predict  a  given  fuel  cell  performance  parameter.  In 
general,  to  allow  for  efficient  use  of  computational  time,  it  is 
customary  in  modeling  to  only  consider  the  details  of  a  given 
transport  process  if  its  characteristic  time  scale  is  within  one 
or  two  orders  of  magnitude  from  the  principle  effect  being 
examined.  Hence,  if  one  is  investigating  load-following  char¬ 
acteristics  of  a  fuel  cell  at  a  time  scale  on  the  order  of  10-2  s, 
then  any  transport  process  having  a  characteristic  time  greater 
than  10°  s  can  be  assumed  constant  over  durations  much  less 
than  10°  s;  hence,  no  detailed  model  equation  is  necessary  to 
relate  the  parameters  of  the  transport  process  to  the  rest  of 
the  model.  On  the  other  hand,  if  there  is  a  transport  process 
having  a  characteristic  time  scale  less  than  10-4  s,  then  the 
related  physical  parameters  can  be  assumed  to  behave  quasi¬ 
steady.  Here,  a  steady  state  equation  is  required  to  relate  the 
parameters  of  the  transport  process  to  the  rest  of  the  model. 

As  an  example  of  time  scale  analysis,  consider  the  case 
of  a  fuel  cell  transient  over  the  duration  of  100  s  with  a  reso¬ 
lution  of  1  s.  Examination  of  the  streamwise  convection  and 
diffusion  in  Table  1,  two  independent  transport  mechanisms 
that  control  streamwise  species  distribution,  shows  that  the 
time  scale  for  convective  transport  is  10_1  s  while  the  time 
scale  for  streamwise  diffusion  is  101  s.  Because  the  convec¬ 
tive  response  is  much  faster  than  the  desired  resolution,  the 
gas  flow  can  be  considered  quasi-steady.  Also,  it  controls  the 
reactant  distribution  much  more  strongly  than  diffusion,  al¬ 
lowing  us  to  ignore  streamwise  diffusion  for  most  analyses 
of  this  given  type  of  fuel  cell  geometry.  Additionally,  the 
time  scale  for  transverse  diffusion  over  the  channel  is  less 
than  10~3  s.  Because  this  is  much  faster  than  the  convec¬ 
tive  time  scale,  the  reactants  in  a  channel  can  be  considered 
largely  uniform  in  the  transverse  direction.  The  ability  to  ig¬ 
nore  these  (and  other)  transport  mechanisms  helps  to  speed 
computational  predictions  making  engineering  analysis  more 
efficient  by  allowing  the  use  of  one-dimensional  quasi-steady 
gas  models,  Fig.  2. 

The  only  time  the  above  method  fails  is  if  the  principle 
effect  being  studied  is  highly  sensitive  to  a  given  transport 
process.  Such  conditions  occur,  for  example,  in  the  analysis 
of  combustion  (e.g.,  the  post-oxidizer  of  a  fuel  cell)  where 
chemical  reactions,  having  time  scales  on  the  order  of  10~3 
to  10-5  s,  occur  (autothermally)  in  a  heated  mixture  of  fuel 
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Fig.  2.  Example  of  planar  co-flow  cell  geometry.  Other  geometries  employ 
similar  to  ID  gas  analysis. 


and  air.  For  these  problems  one  often  sees  that  the  spatial  dis¬ 
tribution  of  the  reaction  (governed  by  convection  time  scale 
of  10°  s)  significantly  depends  on  the  reactions  that  evolve  at 
time  scales  on  the  order  of  10-5  s.  In  fuel  cell  electrochem¬ 
istry,  however,  the  chemical  reactions  are  ‘constrained’  by  the 
external  control  of  the  current  load,  something  that  offers  a 
significant  advantage  over  analysis  in  combustion  chemistry. 

The  present  paper  investigates  two  domains  of  transient 
performance,  and  the  above  characteristic  time  scale  results 
are  used  to  guide  the  development  of  a  proper  model.  The  first 
domain  of  investigation  examines  the  ability  of  the  model  to 
capture  the  fast  transient  current  loading  behavior  of  a  small 
‘button’  fuel  cell  down  to  a  resolution  of  about  0.005  s.  The 
second  domain  of  investigation  applies  the  model  to  a  full 
cell  (10  cm  x  10  cm),  and  is  used  to  predict  the  transient  per¬ 
formance  of  the  fuel  cell  at  a  resolution  of  about  10  s  fol¬ 
lowing  a  load  change.  Hence,  the  former  investigation  will 
include  electrode  diffusion  effects,  while  the  latter  investiga¬ 
tion  will  limit  its  consideration  to  the  thermal  response  and 
treat  the  gas  phase  as  quasi- steady.  Both  cases  will  follow  the 
quasi-steady  electrochemical  response  that  results  from  their 
respective  load  changes. 

2.2.  Model  assumptions 

The  analysis  considers  the  behavior  of  a  single  unit  cell 
located  within  an  infinite  series  of  cells  undergoing  the  same 
dynamics.  Hence,  the  overall  thickness  of  the  unit  cell  (sepa¬ 
rator  plate  to  separator  plate  (also  called  interconnect  plate)) 
is  a  symmetry  length,  and  for  computational  efficiency,  is  the 
only  length  that  needs  to  be  resolved.  Here,  the  term  cell  is 
used  to  describe  the  lumped  fuel  electrode  and  air  electrode 
together  with  the  electrolyte  material. 

Key  assumptions  of  the  model  are: 

1 .  One-dimensional  behavior  along  the  streamwise  direction 
of  all  gas  channels. 

2.  In-plane  heat  conduction  can  be  neglected — an  order-of- 
magnitude  analysis  for  the  Peclet  number  (ratio  of  in¬ 
plane  convective  transport  to  solid  conduction)  shows 
that  convection  in  the  gas  channels  are  4-10  times 
larger  than  conduction,  even  for  composite  nickel  fuel 
electrode  supported  geometries  having  conductivity  of 
~6.2  Wm-1  K"1  [14]. 

3.  Ideal  gas  behavior  with  variable  specific  heats. 

4.  Fumped  temperature  for  the  cell  and  separator  plates  along 
the  transverse  direction  due  to  low  Biot  number  of  the  cell 
convection  (Bi  =  htlk-xflx l  in  Table  1). 

5.  Gas-to-solid  heat  transfer  is  by  convection.  Radiation  is 
ignored  at  this  time,  but  in  general  can  be  important. 

6.  There  is  no  carbon  deposition. 

7.  No  gas  phase  reactions  occur. 

8.  Fuel  electrode  electrochemical  activation  overpotential  is 
negligible  [15]. 

9.  The  fuel  cell  has  100%  current  efficiency — all  reactants 
generate  their  ideal  number  of  electrons. 
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Table  2 


Fundamental  one-dimensional  dynamic  equations  for  an  ideal  gas 


A  dCfc  1  dCfcVx  £>/ 

dt  1  dx  nsk 

(FI) 

*"■  * 

II 

1 

OjI 

*  [*0 

1 

(F2) 

Mw  -  =  gs  +  FL  -  Eb 

(F3a) 

=  QS  +  FL  +  Eb 

(F3b) 

Primarily  due  to  assumption  6  above,  the  model  to  be  pre¬ 
sented  is  most  accurate  for  fuel  electrode  conditions  having 
H2  +  H2O  with  little  or  no  carbon  species.  This  is  acceptable 
for  the  present  work  given  that  the  goal  is  to  provide  a  ‘first- 
look’  at  some  of  the  unsteady  performance  results  that  can 
occur  in  SOFC  systems. 

2.3.  Gas  phase  equations 

The  fundamental  dynamic  equations  required  for  solution 
of  the  one-dimensional  flow  in  a  fuel  cell  are  listed  in  Table  2. 
Here,  the  total  derivative,  D/D t  includes  only  the  convective 
transport  in  the  V  direction,  as  shown  explicitly  for  Eq.  (FI). 
Eqs.  (FI),  (F2),  (F3a)  and  (F3b)  describe  the  conservation  of 
species,  momentum,  and  thermal  energy,  respectively  [16]. 
In  Table  2,  A  is  the  flow  area  of  the  one-dimensional  passage, 
Ck  and  R$k  are  the  mole  density  and  net  surface  mole  supply 
rate  (mol  (v-unit  distance)-1)  of  species  k,  respectively,  vx  the 
velocity  in  the  v-direction,  P  the  pressure,  Mb  the  momentum 
lost  at  the  surface  due  to  flow  at  the  boundaries,  p  the  mass 
density,  Fx  the  net  shear  force  in  the  v-direction  on  the  side 
surfaces  of  the  elemental  control  volume,  Fl  the  thermal  en¬ 
ergy  generated  due  to  fluid  shear  force,  h  the  enthalpy,  E\y  the 
enthalpy  addition  due  to  flow  at  the  side  boundaries,  and  Qs 
the  surface  heat  into  the  elemental  control  volume. 

Multiplying  the  left-  and  right-hand  sides  of  Eq.  (F2)  by  v 
and  adding  (F3a)  and  (F3b)  results  in  an  equation  for  the  con¬ 
servation  of  total  (thermal  and  mechanical)  energy.  Taking 
these  latter  equations  and  Eqs.  (FI)  and  (F2)  and  integrating 
over  the  v-direction  of  a  single  control  volume  (Fig.  2)  results 
in  the  following  set  of  dynamic  equations  to  be  solved  for  each 
finite  control  volume  that  is  present  in  a  one-dimensional  ar¬ 
ray  of  control  volumes  (see  also  [4]).  For  species  conserva¬ 
tion: 

vjt  =  N -  Ni°* +  Rk  (1) 

For  momentum  conservation: 
dpv 

F  —  —  Minlet  Mr  +  Pinlet^I  inlet 

dt 

^exit^exit  T  Fs  (2) 

For  total  energy  conservation: 

dpE  .  . 

F  —  —  ^mlet^inlet  l^exit^exit  T  T  Qs  +  P  (3) 

dt 

Definitions  for  the  various  parameters  are  given  in  nomen¬ 
clature.  The  subscripts  ‘inlet’  and  ‘exit’  denote  the  left  and 


Table  3 


Cell  data  and  operating  conditions 


Parameter 

Unit 

Value 

(button  cell) 

Value 

(planar  cell) 

Fuel  electrode  specification 

Height  of  flow  channel,  h 

mm 

— 

1.0 

Electrode  thickness,  fe 

mm 

1.025 

1.0 

Porosity/tortuosity 

— 

See  text 

0.11 

Temperature,  T 

K 

1073 

1023 

Pressure,  P 

Pa 

101100 

Same 

Supply  H2  mole  fraction 

— 

0.85 

0.54 

Supply  H2O  mole  fraction 

— 

0.15 

0.16 

Supply  CO2  mole  fraction 

— 

— 

0.30 

Electrode  heat  capacity,  Cp 

p-r 

CTQ 

1 

R 

1 

640 

Same 

Electrode  density,  p 

kgm  3 

4200 

Same 

H2  supply,  Sn2 

slpm 

0.050 

1.2 

Air  electrode  specification 

Height  of  flow  channel,  h 

mm 

— 

2.0 

Electrode  thickness,  fe 

[xm 

68 

50 

Porosity/tortuosity 

— 

See  text 

0.11 

Active  area 

2 

cm 

2.19 

100 

Inlet  temperature,  T 

K 

1073 

1023 

Pressure,  P 

Pa 

101100 

Same 

Supply  O2  mole  fraction 

— 

0.204 

0.21 

Supply  N2  mole  fraction 

-0.766 

0.79 

Supply  H2O  mole  fraction 

— 

0.03 

0.0 

Electrode  heat  capacity,  Cp 

Jkg-1  K-1 

570 

Same 

Electrode  density,  p 

kgm-3 

3300 

Same 

Exchange  current  density,  io 

Am-2 

285 

3250 

Transfer  coefficient,  a 

— 

0.5 

0.5 

Air  supply,  SAk 

slpm 

1.0 

5.5 

Electrolyte  specification 

Electrolyte  thickness,  t 

[xm 

7 

5 

Electrolyte  heat  capacity,  Cp 

«— 1 

CTQ 

1 

R 

1 

600 

Same 

Electrolyte  density,  p 

kgm  3 

5100 

Same 

Separator  specification 

Thickness,  t 

mm 

NA 

0.7 

Separator  heat  capacity,  Cp 

Jkg-1  K-1 

NA 

450 

Separator  density,  p 

kg  m~3 

NA 

8900 

Cell  specification 

Contact  resistance,  Rc 

Q  cm2 

0.04 

0.06 

Ohmic  resistance,  R0 

£2  cm2 

0.105 

0.10 

right  surfaces  of  the  control  volume.  The  molar  source  of 
species  k  for  a  given  control  volume,  74,  arises  from  the  elec¬ 
trochemistry  occurring  at  the  surface  of  the  electrolyte.  Under 
normal  fuel  cell  operation  (/>0Acm-2),  the  electrochem¬ 
istry  consumes  O2  on  the  air  electrode  and  H2  on  the  fuel 
electrode ,  and  produces  H2O  on  th z  fuel  electrode}  Under 
reverse  fuel  cell  operation  (electrolysis  mode,  i  <  0  A  cm-2), 
the  electrochemistry  produces  O2  on  the  air  electrode  and  H2 
on  the  fuel  electrode  and  consumes  H2O  on  the  fuel  electrode. 
These  sources  74,  are  directly  related  to  the  electric  current 
according  to  the  number  of  electrons,  %,  associated  with  each 


1  Because  this  paper  examines  conditions  of  both  normal  and  reverse  fuel 

cell  operation  modes,  the  terms  ‘cathode’  and  ‘anode’  become  confused, 

since  depending  on  the  mode  of  operation  they  will  apply  to  opposite  elec¬ 
trodes.  Instead,  the  paper  will  use  the  terms  ‘air  electrode’  and  ‘fuel  elec¬ 
trode’  to  represent  the  electrodes  where  O2  and  H2  +  H2O  exist,  respectively. 
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specie  molecule  consumed/generated  in  the  electrochemical 
process.  For  the  present  analysis  on  solid  oxide  cells,  four 
electrons  are  processed  by  each  oxygen  molecule  consumed 
at  the  air  electrode,  and  two  electrons  are  processed  by  each 
hydrogen  (and  water  vapor)  molecule  processed  at  the  fuel 
electrode.  Explicitly,  we  have 


rik  for  O2,  H2,  and  H2O,  are  —4,  —2,  and  2,  respectively. 
Hence,  for  i  <  0  A  cm-2,  Ro2  >0  kg  mol  s_1 . 

In  Eq.  (2),  the  momentum  flow,  Mr,  comes  from  the  term 
Mb  in  (F2)  and  is  the  net  momentum  loss  resulting  from  reac¬ 
tant  removal  out  of  the  control  volume  at  the  electrolyte  sur¬ 
face  due  to  the  electrochemical  reaction.  For  mass  addition 
at  the  electrolyte  surface,  it  is  assumed  that  no  net  addition  of 
momentum  occurs — that  is,  there  is  no  favored  direction  to 
the  emission  of  molecules  from  the  electrolyte  due  to  electro¬ 
chemical  reactions,  and  hence  no  net  momentum  is  produced. 

The  parameter,  Lr,  in  Eq.  (3)  is  an  energy  flow  term  that 
resulted  from  Lb  in  Eqs.  (F3).  Lr  represents  the  net  enthalpy 
added  to  the  control  volume  resulting  from  species  emis¬ 
sion/absorption  at  the  electrolyte  surface  due  to  the  electro¬ 
chemical  reaction.  Here,  the  energy  of  the  species  removed 
from  the  gas  phase  is  assessed  at  the  temperature  of  the  gas, 
and  the  energy  of  the  species  added  to  the  gas  phase  is  as¬ 
sessed  at  the  temperature  of  the  cell. 

For  all  control  volumes  except  at  the  exit,  E  is  the  sum  of 
internal  energy  and  kinetic  energy.  For  the  exit  control  volume 
we  assume  pressure  fluctuations  are  zero  (constant  pressure 
boundary  condition),  and  the  enthalpy  plus  kinetic  energy  is 
used.  The  assumption  here  is  that  the  exit  of  the  fuel  cell  goes 
to  near  atmospheric  conditions  (or  some  other  specified  fixed 
exit  pressure)  which  results  in  a  constant  pressure  boundary 
condition.  On  the  other  hand,  the  inlet  pressure  can  vary  with 
time  as  may  be  caused  by  the  action  of  an  upstream  flow 
controller. 

The  model  allows  for  two  types  of  gas  phase  solutions.  One 
is  fully  dynamic  in  that  it  solves  the  above  one-dimensional 
dynamic  equations  for  species,  momentum  and  energy.  The 
other  assumes  quasi- steady  gas  phase  conditions  to  exist  and 
solves  the  steady-state  gas  phase  equations  for  these  param¬ 
eters  given  the  instantaneous  conditions  at  its  boundaries. 
Here,  the  same  equations  listed  above  are  used,  but  the  time- 
derivative  terms  are  dropped  and  the  resulting  algebraic  equa¬ 
tions  are  solved  simultaneously.  Since  flow  times  for  gases 
through  a  cell  are  on  the  order  of  1  s  or  less  (see  Table  1), 
the  quasi-steady  solution  is  useful  to  speed  up  calculations 
on  long  transients  (e.g.,  >10  s). 


2.4.  Solid  phase  energy  equation 

For  control  volumes  present  in  the  cell  and  separator  plate, 
an  energy  equation  similar  to  Eq.  (3)  arises,  but  the  convective 
terms  will  obviously  be  zero.  Also,  for  any  solid,  there  is  an 
additional  energy  source  term  due  to  internal  ohmic  heating. 


For  the  cell  alone,  there  is  a  final  additional  energy  gener¬ 
ation  term  due  to  the  overall  reaction  of  H2  with  O2.  This 
energy  generation  results  from  the  reaction  process  which 
experiences  a  net  change  in  entropy  [17]: 


Q 


rxn  — 


T  AS  AJ 
nF 


where  AS  is  the  change  in  entropy  (per  mole  of  H2)  of  the 
overall  reaction  of  H2  with  O2,  Aa/  the  current  through  the 
control  volume,  n  (no  subscript)  the  number  of  electrons  pro¬ 
cessed  per  mole  of  H2,  and  F  the  Faraday  constant.  Hence, 
the  following  equations  apply  to  the  cell  and  separator  plate, 
respectively: 


d  pE 

V  +  Qs  +  Qrxn  +  iR 

at 


V  -r —  =  Qs  +  iR  (7) 

at 

where  E  is  the  internal  energy  of  the  given  solid  material, 
Lr  the  net  enthalpy  added  to  the  cell  resulting  from  species 
emission/absorption  for  both  fuel  electrode  and  air  electrode 
surfaces,  Qs  the  same  as  before,  and  R  the  cell  net  elec¬ 
tronic  +  ionic  +  contact  cell  resistance  (Q  m2). 


2.5.  Cell  voltage  and  loss  mechanisms 

The  cell  voltage  is  modeled  by  accounting  for  various  loss 
mechanisms  in  an  otherwise  ideal  cell.  The  ideal  fuel  cell 
voltage  is  given  by 

AG 

En  = - -  (8) 

nF 

AG  is  the  change  in  Gibbs  free  energy  for  reaction  of  H2  with 
O2  at  the  temperature  and  pressure  of  the  cell,  which  provides 
the  maximum  available  work  energy  for  this  reaction.  If  this 
work  energy  is  given  to  n  number  of  electrons,  then  their 
voltage  would  be  L 

Voltage  losses  (overpotentials)  within  a  fuel  cell  result 
from  three  primary  effects:  (1)  electrical  and  ionic  resis¬ 
tance,  Lr;  (2)  change  in  electrolyte  surface  reactant  concen¬ 
trations  from  their  free- stream  values  due  to  diffusion  trans¬ 
port  through  the  electrodes,  Lq;  and  (3)  an  activation  over¬ 
potential  needed  to  drive  the  electrochemical  reactions,  La 
[18].  As  shown  below,  these  losses  depend  on  current  flow 
through  the  cell.  The  cell  voltage  under  current  loading  can 
be  defined  as 


Lcell  =  EN  -  Lr  -  ^2  LCi  (9) 

This  cell  voltage  is  the  voltage  that  would  exist  on  some 
external  load  that  draws  current  from  the  cell. 

The  net  cell  electronic  +  ionic  +  contact  resistive  losses  are 
modeled  as 


Lr  =  Ri  (10) 

where  i  is  the  local  current  density  for  a  given  control  volume. 
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For  diffusive  losses,  the  model  can  either  employ  a  bulk 
quasi-steady  diffusion  loss  mechanism,  or  resolve  the  dy¬ 
namic  diffusion  through  the  electrodes  for  each  control  vol¬ 
ume.  The  latter  model  is  given  in  detail  by  Gemmen  et  al. 
[19]. 

The  quasi-steady  diffusion  loss  model  uses  the  following 
set  of  equations  derived  from  an  analysis  of  the  effect  of  diffu¬ 
sion  on  reactant  concentrations  near  the  electrolyte/electrode 
interface  and  the  subsequent  loss  to  the  cell  Nernst  potential. 
For  reactant  k  diffusion  toward  the  electrolyte  interface,  the 
loss  is  given  as 

RUT  /1-A 

La  = - —  In  (  — —  )  (1  la) 

nF  \  iLk  ) 


where  i\±  is  the  diffusion  limiting  current  density  for  reactant 
k ,  which  is  modeled  as 


ihk 


CkoDeffnF 

L 


(12a) 


where  Deff  is  the  effective  diffusivity  for  the  transported  re¬ 
actant,  and  tQ  the  electrode  diffusion  thickness.  (Note  that 
because  of  the  formulation  given  above,  the  sign  of  z'l  is 
always  the  same  as  i.  This  understanding  is  important  when 
later  we  discuss  conditions  for  current  reversal  in  a  cell,  i  <  0.) 
For  reactant  diffusion  away  from  the  electrolyte  interface,  the 
quasi-steady  loss  model  is 


RUT 


nF 


In  1  - 


(Hb) 


where  now  the  limiting  current  density  is  given  by 


ihk 


Q)(l  yko)F>effnF 

1q 


(12b) 


circuit  at  each  electrode,  where  the  resistor  represents 
the  electrochemical  loss  mechanism.  The  model  can  be 
optionally  set  to  solve  for  the  double-layer  transient,  or 
to  assume  fully  quasi-steady  loss  behavior  (Eq.  (14)).  To 
solve  for  the  transient  double-layer  capacitance  behavior, 
the  model  employs  the  following  equation  [21]: 


d?)  _  ic 

d  t  Ci  i 


(15) 


where  ic  is  the  current  density  that  charges  the  double-layer 
capacitance  of  Cdi  (F  m~2).  A  generally  recognized  value  for 
Cdi  has  not  been  established,  but  values  as  high  as  1  F  m-2 
have  been  suggested  [21 ,22] .  Because  of  the  parallel  R-C  cir¬ 
cuit,  the  instantaneous  current  through  a  given  computational 
node  is  the  sum  of  current  through  both  R  and  C  elements. 
Hence, 


ic 


i  —  l  R 


(16) 


with  the  assumption  of  quasi-steady  electrochemical  loss, 
z'r  is  given  by  Eq.  (14).  Because  the  time  scale  for  the 
double-layer  capacitance  is  on  the  order  of  100  p,s  is  or 
less  (see  Table  1),  the  present  work  uses  the  quasi-steady 
Butler- Volmer  model. 

Note  that  because  the  fuel  and  oxidant  are  consumed  as 
they  pass  through  the  cell,  the  ideal  cell  voltage  at  the  exit 
will  be  lower  than  that  at  the  inlet.  The  actual  overall  cell 
voltage  achieved  will  approach  the  lowest  achievable  voltage 
due  to  the  high  conductivity  of  the  electrode  surfaces  [1 8] .  As 
usual  for  loss  mechanisms,  it  is  assumed  that  the  unachieved 
voltage  over  the  surface  is  dissipated  as  heat  directly  into  the 
cell  material. 


2.6.  Reversible  fuel  cell  operation 


For  positive  current  (z>0.0),  Eqs.  (11a)  and  (12a)  apply  to 
O2  on  the  air  electrode  and  H2  on  the  fuel  electrode,  while 
Eqs.  (lib)  and  (12b)  apply  to  H2O  on  the  fuel  electrode.  For 
negative  current,  Eqs.  (lib)  and  (12b)  apply  to  O2  on  the  air 
electrode  and  H2  on  the  fuel  electrode,  while  Eqs.  (11a)  and 
(12a)  apply  to  H2O  on  the  fuel  electrode. 

For  the  activation  loss,  La,  we  have 


La  =  0c  +  ?7a 


(13) 


As  shown  in  Eq.  (13),  in  general  there  is  an  activation  loss  at 
both  air  and  fuel  electrodes,  ijc  and  77a,  respectively.  For  the 
work  here,  because  air  electrode  losses  dominate  La,  the  fuel 
electrode  activation  is  ignored. 

Under  steady  state,  the  electrochemical  loss  is  often  rep¬ 
resented  using  the  Butler- Volmer  equation  [20] : 


z  =  z0 


—  exp 


^(1  —  a)  nFrj^j 


(14) 


Because  of  a  double-layer  capacitance  at  the 
electrolyte-electrode  interface,  there  is  actually  a  short 
transient  behavior  to  the  activation  loss.  The  transient  loss 
is  typically  represented  as  a  parallel  resistor-capacitance 


The  models  for  the  loss  terms  cited  above  are  common  to 
fuel  cell  modeling  [8].  However,  most  past  modeling  efforts 
simply  consider  conditions  where  current  flow,  z,  through  the 
cell  is  positive.  As  will  be  shown  below,  there  are  transient 
conditions  where  reverse  current,  z  <  0,  may  arise.  These  con¬ 
ditions  arise  when  Lcen  is  greater  than  Ln,  which  means  that 
the  loss  terms  on  the  right-hand  side  of  Eq.  (9)  must  become 
positive.  (Clearly,  to  force  reverse  current,  the  operating  elec¬ 
trical  potential  must  be  greater  than  the  chemical  (i.e.,  Nernst) 
potential.)  Because  rate  limiting  processes  are  present  under 
these  conditions  just  as  for  positive  current  conditions,  poten¬ 
tial  losses  occur,  and  using  the  above  models  requires  these 
losses  to  be  positive- valued.  Given  the  formulations  above, 
proper  numerical  treatment  for  both  resistance  and  activation 
losses  (Eqs.  (10)  and  (14))  occurs  automatically  for  both  pos¬ 
itive  and  negative  current  domains,  but  for  proper  accounting 
of  diffusion  under  negative  current  conditions,  the  sign  on 
diffusion  terms  (Lq)  must  be  reversed.  This  can  be  accom¬ 
modated  by  changing  Eq.  (9)  to  be 

Lceii  =  En  ~  Lr  +  ^2  La  ~  Fa  (17) 

for  negative  current  conditions. 
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3.  Numeric  representation  and  solution  technique 

The  above  equations  allow  the  transient  solution  of  one¬ 
dimensional  arrays  of  both  gas  and  solid  bodies  that  comprise 
a  fuel  cell.  Depending  on  the  needs  of  the  engineer,  different 
submodels  can  be  employed  (e.g.,  resolving  the  double-layer 
transient  may  not  be  needed  for  studies  looking  at  the  longer 
time  scales  associated  with  cell  thermal  transients  and  their 
impact  on  cell  voltage).  Finite  control  volumes  for  both  gas 
passages  and  solid  bodies  are  identified  and  properly  associ¬ 
ated  to  provide  proper  conservation  of  mass,  momentum  and 
energy.  The  NETL  code  that  implements  the  above  equations 
can  solve  for  co-flow,  counter-flow,  and  cross-flow  fuel  cell 
geometries.  For  the  cross-flow  geometry,  an  array  of  fuel  elec¬ 
trode  channels  and  air  electrode  channels  are  generated  and 
assembled  to  form  such  a  cell,  similar  to  the  work  of  Ahmed 
et  al.  [13].  For  the  co-flow  and  counter-flow  geometries,  only 
one  channel  of  both  fuel  electrode  and  air  electrode  gas  is 
required. 

Any  number  of  control  volumes  are  possible  along  the 
streamwise  direction  for  each  of  the  four  major  components: 
fuel  electrode  gas,  air  electrode  gas,  cell,  and  separator  plate. 
For  the  work  reported  here,  eight  control  volume  ‘nodes’  are 
used  which  provides  sufficient  resolution  to  capture  fuel  cell 
performance  [12].  For  solving  dynamic  gas  phase  problems, 
the  equations  are  solved  following  the  techniques  outlined  by 
Patankar  [23],  wherein  the  pressure  and  velocity  nodes  are 
staggered. 

In  the  solution  of  the  model,  either  the  cell  voltage  or  total 
cell  current  is  specified  by  the  user  (voltage  control  versus 
current  control).  To  solve  conditions  of  the  cell  at  every  time 
step  based  on  specified  instantaneous  boundary  conditions, 
the  model  first  solves  for  all  node  electrochemical  reaction 
rates  (current  densities)  given  the  specified  cell  voltage  or  to¬ 
tal  current.  If  a  quasi-steady  activation  loss  is  assumed  (i.e., 
no  double-layer  dynamics  to  be  solved),  then  an  iterative  ap¬ 
proach  is  used  to  determine  the  cell  current  density  profile. 
Node  current  at  each  time  step  is  iterated  so  as  to  ensure  a 
uniform  cell  voltage  (to  within  4  piV).  A  uniform  cell  voltage 
is  a  common  assumption  in  fuel  cell  modeling  given  the  high 
electric  conductivity  of  the  electrodes  and  separator  plates.  A 
combination  between  a  Newton  and  simple  bisection  method 
is  used  to  converge  to  the  desired  solution.  Once  the  current 
is  known,  the  dynamic  solution  for  mass,  momentum  and  en¬ 
ergy  (Eqs.  (l)-(3))  for  all  cell  domains  is  determined  through 
an  explicit  first-order  backward  Euler  method.  For  a  quasi¬ 
steady  gas  solution,  the  conservation  equations  for  mass,  mo¬ 
mentum,  and  energy  are  solved  simultaneously  in  a  coupled 
fashion  with  the  electrochemical  solution. 


4.  Model  application 

The  model  has  been  tested  for  self-consistency  regarding 
the  application  of  its  various  submodels.  For  example,  the 
quasi- steady  gas  solution  was  tested  and  found  to  provide  the 
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Fig.  3.  Experimental  hardware. 

same  conditions  as  the  fully  dynamic  gas  model  predictions 
when  taken  to  steady  state.  The  dynamic  model  was  also 
verified  by  comparing  its  results  to  that  of  the  National  Fuel 
Cell  Research  Center  and  found  to  be  comparable  [24] . 

The  model  was  also  tested  in  cross-flow  geometry  by  com¬ 
paring  results  to  that  of  Ahmed  et  al.  [13] — the  results  were 
the  same  to  within  4%.  In  the  present  work,  we  report  the 
success  of  model  validation  using  experimental  cell  transient 
data.  The  cell  transient  data  was  obtained  from  a  button  cell 
operated  inside  an  enclosed  rig  that  controlled  both  the  flow 
of  air  and  fuel.  A  view  of  the  test  set  up  is  shown  in  Fig.  3. 

Finally,  the  model  was  applied  to  investigate  three  cases  of 
planar  10  cm  x  10  cm  fuel  cell  design  (CoF,  CrF  and  CtF).  For 
each  case,  several  load  transients  were  imposed  that  are  repre¬ 
sentative  of  real  world  applications  where  a  system  might  ex¬ 
perience  anything  from  small  changes  in  load  (<30%)  to  com¬ 
plete  load  loss.  For  simplicity,  following  the  work  of  Ahmed 
et  al.  [13],  eight  channels  for  fuel  electrode  flow  and  eight 
channels  for  air  electrode  flow  are  considered  using  the  same 
one-dimensional  model  described  above  for  each  channel. 

4.1.  Button  cell  study 

For  the  transient  button  cell  cases,  a  fully  dynamic  gas 
solution  was  calculated  using  the  electrode  diffusion  model 
of  Gemmen  et  al.  [19],  but  the  double-layer  transient  was 
ignored  (which  again  has  a  time  constant  on  the  order  of 
100  pus  and  is  very  short  relative  to  the  transient  under  study). 
Given  the  uniformity  of  the  conditions  over  the  cell  that  can 
be  expected  from  the  known  flow  conditions  used  in  the  ex¬ 
periment,  the  model  uses  a  single  control  volume  for  the  fuel 
electrode  and  air  electrode  systems.  The  region  up  to  3  cm 
from  the  surface  of  the  cell  was  included  in  the  model  anal¬ 
ysis.  Load  perturbation  occurred  through  current  control. 

Values  for  the  model  parameters  used  in  the  button  cell 
comparisons  are  shown  in  Table  1 .  The  values  for  the  param¬ 
eters  were  determined  from  the  cell  specifications  given  by 
the  University  of  Utah  from  where  the  cells  were  obtained, 
and  from  the  known  conditions  of  the  test  fixture.  The  actual 
button  cell  included  thin  fuel  electrode  and  air  electrode  in¬ 
terlayers,  25  and  18  gm,  respectively,  and  these  thicknesses 
were  incorporated  into  the  bulk  electrode  to  establish  a  single 
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Fig.  4.  Model  node  ID  numbers  used  on  charts. 


electrode  for  the  model.  This  assumption  should  be  reason¬ 
able  given  the  relatively  thin  nature  of  these  layers.  The  con¬ 
tact  resistance  value  shown  (obtained  by  matching  the  exper¬ 
imental  VI  data)  was  attributed  to  the  experimental  method 
used  to  contact  the  fuel  electrode  and  air  electrode  current 
leads  to  the  button  cell. 

4.2.  Planar  cell  cases 

For  the  three  CoF,  CrF  and  CtF  transient  studies,  the  quasi¬ 
steady  gas  and  activation  loss  models  were  used  together  with 
the  bulk  diffusion  model,  and  the  model  was  used  to  predict 
the  response  of  the  cell  to  load  transients.  This  study  em¬ 
ployed  voltage  controlled  transients.  The  model  parameters 
used  for  these  studies  are  also  shown  in  Table  3.  These  param¬ 
eter  values  provide  cell  performance  close  to  what  is  antici¬ 
pated  for  SECA  Program  fuel  cell  systems — at  a  cell  voltage 
of  0.77  V  and  fuel  and  air  utilizations  of  0.81  and  0.11,  re¬ 
spectively,  cell  average  current  density  ~0.74Acm~2.  The 
cells  were  assumed  to  have  a  10  cm  x  10  cm  active  area,  so 
the  channel  lengths  for  both  fuel  electrode  and  air  electrode 
are  10  cm  long.  As  mentioned  previously,  eight  channels  are 
used  for  both  the  fuel  electrode  and  air  electrode,  and  each 
channel  has  eight  nodes.  The  resultant  size  of  each  node  is 
12.5  mm  x  12.5  mm.  For  reference,  the  node  (control  vol¬ 
ume)  ID  numbers  for  the  model  are  shown  in  Fig.  4.  Time 
steps  used  in  the  solution  of  the  model  ranged  from  3.0E— 2  to 
8.0E— 2  s,  which  were  controlled  by  the  solution  of  the  solid 
phase  thermal  transient  (the  gas  was  assumed  quasi- steady). 


5.  Results 

Results  for  the  button  cell  and  10  cm  xlOcm  cell  cases 
are  now  presented. 

5.1.  Button  cell  case 

Transient  button  cell  experimental  and  model  results  are 
shown  in  Fig.  5.  Two  load  changes  were  examined,  one  with 
load  reduced  from  1  to  0A,  the  other  from  2  to  0A.  For 
each  case,  the  data  was  taken  at  both  a  short  time  (fast  data 
speed)  and  a  long  time  (slow  data  speed)  in  order  to  capture 


Fig.  5.  Transient  button  cell  data  and  model  comparison.  Inset  shows  same 
data  over  shorter  time  scale. 

the  transient  as  closely  and  fully  as  possible.  As  can  be  seen, 
the  two  sections  of  the  transient  data  align  very  well,  indicat¬ 
ing  good  repeatability  of  the  experimental  data.  Results  for 
the  1  A  load  change  employed  a  porosity /tortuosity  value  of 
0.025,  while  the  2  A  load  change  employed  a  value  of  0.020. 
These  were  selected  to  best  fit  the  data.  The  model  results 
agree  fairly  well  with  the  experimental  data,  which  provides 
some  measure  of  confidence  in  the  model. 

5.2.  Cross -flow  case — load  increase 

Results  from  the  cross-flow  study  examining  load  increase 
dynamics  are  shown  in  Figs.  6  and  7.  Fig.  6  shows  transient 
current  densities  at  several  node  locations  on  the  cell  for  a 
0.96-0.84  V  load  change  case — see  Fig.  4  for  a  key  to  the 
node  ID  numbers.  As  can  be  seen,  there  is  an  initially  high 
current  flow  through  the  cell,  and  then  as  the  temperature  of 
the  cell  slowly  increases  due  to  internal  heating,  the  Nernst 
voltage  gradually  drops  which  causes  a  slow  drop  in  current 
flow.  It  is  also  seen  that  while  the  cathode  inlet  and  exit  nodes 
show  fairly  similar  current  density  before  and  initially  after 
the  load  change,  once  the  cell  has  reached  its  new  equilibrium 
following  the  load  change,  the  inlet  and  exit  nodes  are  signif¬ 
icantly  different.  Finally,  it  is  seen  that  an  overall  increase  in 
current  by  a  factor  of  about  3.75  results  from  this  load  change. 
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Fig.  6.  Cross-flow  load  increase  (0.96-0.84  V). 
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Fig.  7.  Cross-flow  load  increase  (0.96-0.70  V). 


Fig.  9.  Cross-flow  load  decrease  (0.70-0.96  V). 


Fig.  7  shows  similar  results  as  in  Fig.  6  but  for  a  larger 
0.96-0.70  V  load  change.  For  this  transient,  there  is  an  overall 
current  increase  by  a  factor  of  about  6.6.  While  most  of  the 
results  highlighted  for  the  prior  case  of  Fig.  6  are  present  again 
here,  this  case  shows  a  distinctly  different  transient  result  for 
the  air  electrode  exit  nodes  (e.g.,  57  and  64).  This  difference 
will  be  discussed  in  more  detail  below. 

5.3.  Cross -flow  case — load  decrease 


new  equilibrium  state.  Here,  an  overall  current  decrease  by  a 
factor  of  about  0.87  results. 

Finally,  Fig.  10  shows  a  contour  plot  of  the  current  density 
for  this  same  case  immediately  following  the  load  change. 
As  can  be  seen,  a  significant  portion  of  the  cell  experiences 
this  negative  current  state  following  the  strong  load  decrease. 
Close  examination  of  the  model  data  at  this  time  show  that 
the  peak  negative  current  density  achieved  on  the  cell  is  about 
—  1330  Am-2.  The  average  cell  load  current  at  this  time  is 
about  —9  Am-2. 


Results  from  the  cross-flow  study  examining  load  decrease 
dynamics  are  shown  in  Figs.  8  and  9.  Fig.  8  again  shows 
transient  current  densities  at  several  node  locations  on  the 
cell  but  now  for  a  0.70-0.77  V  load  change.  The  results  show 
a  reversed  transient  as  compared  to  the  load  increase  case  as 
might  be  expected.  In  particular,  there  is  a  sudden  decrease 
in  current  flow,  which  is  followed  by  a  gradual  increase  in 
current  as  the  temperature  of  the  cell  slowly  drops  allowing 
the  Nernst  potential  to  increase  which  drives  more  current. 
For  this  case,  there  is  an  overall  current  decrease  by  a  factor 
of  about  0.25. 

Fig.  9  shows  similar  results  but  for  a  0.7-0.957  V  load 
change.  These  results  show  an  additional  feature  not  seen  in 
the  prior  results,  namely  negative  current  for  several  nodes 
over  an  extended  period  of  time  until  the  cell  has  reached  a 
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Fig.  8.  Cross-flow  load  decrease  (0.7-0.77  V). 
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5.4.  Co-flow  case — load  decrease 

Similar  load  change  cases  were  studied  for  the  co-flow 
geometry  as  for  the  above  counter-flow  geometry.  Because 
of  the  unique  results  shown  under  strong  load  decreases,  we 
present  only  results  for  the  0.7-0.957  V  load  change.  The 
results  for  the  other  load  change  cases  in  general  followed 
the  cross-flow  geometry  results.  Fig.  1 1  shows  a  contour  plot 
of  the  current  density  immediately  following  the  load  change. 
Again,  a  significant  portion  of  the  cell  experiences  negative 
current  conditions  over  a  lengthy  duration  as  the  cell  slowly 
thermally  equilibrates  to  its  new  steady  state  condition.  Close 
examination  of  the  model  data  at  this  time  show  that  the 
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Fig.  1 1 .  Contour  showing  region  of  current  reversal  following  load  decrease 
(0.70-0.96  V)  (co-flow). 

peak  negative  current  density  achieved  on  the  cell  is  about 
—690 Am-2.  The  average  cell  load  current  at  this  time  is 
about  —74  Am-2. 

5.5.  Counter-flow  case — load  decrease 


6.1.  Button  cell  case 

The  close  match  between  model  results  and  experimen¬ 
tal  data  shown  in  Fig.  5  are  encouraging.  The  need  for  the 
slightly  different  tortuosity/porosity  values  between  the  two 
cases  (in  order  to  best  match  the  data)  need  to  be  investi¬ 
gated  further.  Nonetheless,  these  early  results  suggest  that 
the  model  properly  accounts  for  most  of  the  necessary  con¬ 
trolling  physical  issues.  It  also  suggests  that  the  method  em¬ 
ployed  here  could  be  used  to  help  assess  some  of  the  important 
transport  parameters  of  cells.  For  example,  one  can  perform 
such  high  speed  transient  measurements  and  then  use  these 
models  to  determine  correct  values  of  the  parameters  that 
control  diffusion,  such  as  effective  tortuosity/porosity  which 
are  often  difficult  to  assess,  especially  for  thin  layers.  How¬ 
ever,  more  work  needs  to  be  done  to  determine  general  va¬ 
lidity  or  applicability.  In  particular,  experimental  data  on  the 
transient  performance  of  full  10  cm  x  10  cm  cells  would  be 
helpful. 

6.2.  Cross -flow  case 


Similar  load  change  cases  were  again  studied  for 
the  counter-flow  geometry.  Again,  only  results  for  the 
0.7-0.957  V  load  change  are  presented.  The  results  for  the 
other  load  change  cases  in  general  followed  the  cross-flow 
geometry  results.  Fig.  12  shows  a  contour  plot  of  the  cur¬ 
rent  density  immediately  following  the  load  change.  Again, 
a  significant  portion  of  the  cell  experiences  negative  current 
conditions  over  a  lengthy  duration  as  the  cell  slowly  thermally 
equilibrates  to  its  new  steady  state  condition.  Close  examina¬ 
tion  of  the  model  data  at  this  time  show  that  the  peak  negative 
current  density  achieved  on  the  cell  is  about  —1060  Am-2. 
The  average  cell  load  current  density  at  this  time  is  about 
350  Am-2. 


6.  Discussion 

The  results  of  the  button  cell  and  large  cell  studies  are 
discussed  in  further  detail  in  the  following  sections. 
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Fig.  12.  Contour  showing  region  of  current  reversal  following  load  decrease 
(0.70-0.96  V)  (counter-flow). 


Application  of  the  model  to  the  10  cm  x  10  cm  cross-flow 
geometry  shows  distinctly  different  behavior  can  result  de¬ 
pending  on  the  magnitude  and  direction  of  a  load  change. 
For  a  small  load  increase  (Fig.  6),  all  nodes  have  the  same 
characteristic  transient  whereby  there  is  a  sudden  increase  in 
current  followed  by  a  gradual  decay  as  the  cell  equilibrates 
to  its  new  temperature.  On  the  other  hand,  for  a  large  load 
increase  (Fig.  7),  while  most  of  cell  region  toward  the  fuel 
inlet  (e.g.,  nodes  1,8,  17  and  24)  has  a  sharp  increase  in  cur¬ 
rent  immediately  following  a  load  change,  the  region  toward 
the  fuel  exit  (e.g.,  nodes  57  and  64)  has  a  more  gradual  in¬ 
crease  to  its  peak  value.  It  is  also  possible  to  find  a  certain 
point  in  the  cell  where  the  controlling  effects  appear  to  can¬ 
cel  such  that  no  change  in  current  will  occur  (e.g.,  node  49  in 
Fig.  7). 

For  a  small  load  decrease  (Fig.  8),  current  at  all  nodes  show 
largely  the  same  normalized  response.  However,  for  large 
load  decrease  (Fig.  9),  as  may  occur  in  an  emergency-stop 
command,  certain  portions  of  the  cell  exhibit  current  reversal 
as  seen  in  Fig.  10.  The  reversed  current  condition  lasts  for 
approximately  330  s  for  the  case  studied.  It  is  clear  that  over 
this  time  a  certain  portion  of  the  cell  near  the  air  inlet  pro¬ 
vides  enough  current  to  meet  the  need  of  the  external  load  as 
well  as  to  drive  current  in  reverse  over  its  exit  portion.  Further 
analysis  of  the  model  data  shows  that  the  main  driver  for  this 
reverse  current  is  the  cell  thermal  non-uniformity.  It  is  found 
that  the  air  electrode  exit  portion  of  the  cell  results  in  a  lower 
Nernst  potential  as  compared  to  that  at  the  inlet.  Fig.  1 3  shows 
the  Nernst  potential  distribution  on  the  cell  immediately  fol¬ 
lowing  the  load  decrease.  The  decreased  Nernst  voltage  at  the 
exit  is  due  to  the  higher  exit  temperature  (~  1 00  K  higher),  and 
this  is  enough  to  drive  current  in  reverse.  Based  on  the  first- 
order  analysis  for  in  plane  thermal  transport,  improvements  to 
the  model  to  include  such  transport  would  only  slightly  lessen 
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Fig.  13.  Contour  of  Nernst  potential  following  load  decrease  (0.70-0.96  V) 
(cross-flow). 

these  effects  (in  both  time  and  magnitude),  but  will  not  re¬ 
move  them  altogether  since  the  process  is  thermodynamically 
driven. 

To  illustrate  what  is  occurring  during  current  reversal, 
Fig.  14  shows  the  path  of  electric  charge  through  the  cell. 
The  left  portion  of  the  cell  in  the  figure  shows  positive  cur¬ 
rent  conditions  which  consume  O2  and  H2  and  produce  H2O 
at  the  fuel  electrode  surface.  The  transferred  electrons  are 
provided  to  the  right  side  of  the  fuel  electrode  where  H2O  is 
reduced  to  generate  O2-  ions  which  travel  in  reverse  through 
the  electrolyte  to  generate  O2  at  the  air  electrode.  The  elec¬ 
trons  released  on  the  right  side  of  the  air  electrode  traverse 
back  to  the  left  portion  of  the  cell  and  the  process  con¬ 
tinues  until  the  thermodynamic  driving  potential  no  longer 
exists  to  drive  this  process.  As  a  check  for  consistency  in 
the  numeric  calculation,  Figs.  15  and  16  show  the  hydrogen 
concentration  that  results  just  prior  to  and  after  the  unload 
event,  respectively.  On  the  whole,  Fig.  15  shows  a  gradual 
decay  in  hydrogen  concentration  as  the  fuel  electrode  flow 
passes  through  the  cell,  as  expected.  Fig.  16  shows  a  near 
recovery  of  the  hydrogen  concentration  to  values  close  to 
the  supply  concentration  (yH2  =  0.54),  as  would  be  expected 
for  an  unload  event.  However,  they  are  modified  slightly  by 
the  consumption  of  hydrogen  toward  the  air  electrode  in¬ 


Fig.  15.  Hydrogen  concentration  prior  to  unload  event. 

let,  and  the  creation  of  hydrogen  at  the  air  electrode  exit 
(per  the  current  density  field  shown  in  Fig.  10).  The  self- 
consistency  of  these  results  confirms  that  the  current  rever¬ 
sal  identified  is  not  likely  due  to  some  anomalous  numerical 
computation. 

6.3.  Co-flow ;  counter-flow  and  cross -flow  comparison 

The  temperature  profile  has  been  identified  as  the  gov¬ 
erning  driving  method  for  the  reverse  current  following  load 
decrease.  As  such,  the  CoF  and  CtF  results  can  be  expected 
to  be  different  from  the  CrF  results,  and  follow  the  thermal 
characteristics  of  those  designs.  This  result  is  clearly  evi¬ 
dent  in  the  current  density  contours  shown  in  Figs.  1 1  and  12 
for  the  CoF  and  CtF  cases,  respectively.  Past  work  by  oth¬ 
ers,  Khaleel  et  al.  [25],  have  already  shown  that  CtF,  CrF,  and 
CoF  geometries  result  in  different  temperature  variations  over 
the  cell,  d77T=  0.352,  0.349,  0.241,  respectively.  The  results 
here  for  current  reversal  are  consistent  with  these  thermal  re¬ 
sults  whereby  the  CtF,  CrF  and  CoF  geometries  show  a  cur¬ 
rent  reversal  from  the  average  current  (R  =  Min(/)  —  Ave(/)) 
as  —1410,  —1321,  —616  Am-2,  respectively.  It  is  apparent 
that  the  overall  temperature  difference  is  a  sufficient  indica¬ 
tor  for  predicting  the  presence  of  current  reversal  during  high 
unload  events. 
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Fig.  14.  Schematic  representation  of  electron  charge  path  occurring  during 
reverse  current. 


Fig.  16.  Hydrogen  concentration  after  unload  event. 
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The  aforementioned  cases  all  employed  a  0.7-0.957  V 
load  change.  The  model  was  implemented  a  second  time  for 
each  geometry  to  determine  the  acceptable  limit  of  load  varia¬ 
tion  (in  terms  of  current)  that  would  avoid  the  current  reversal 
conditions.  For  the  CrF,  CoF  and  CtF  geometries,  the  limit  of 
load  current  decrease  was  83,  85  and  86%,  respectively,  from 
their  full  load  cell  voltage  condition  of  0.7  V.  Load  changes 
that  would  induce  a  greater  decrease  in  cell  load  current  will 
cause  at  least  some  portion  of  the  cell  to  produce  reversed 
current. 

The  reverse  current  conditions  identified  by  the  model 
have  not  been  described  previously,  and  more  work  is  needed 
to  examine  the  limiting  conditions  where  such  reversal  oc¬ 
curs.  If  such  conditions  are  indeed  found  to  be  possible  for  a 
given  fuel  cell  system,  then  the  SOFC  system  designer  may 
need  to  consider  its  impact  on  the  cell  components.  For  ex¬ 
ample,  it  is  not  yet  clear  if  the  thermodynamic  potential  that 
drives  this  current  reversal  will  result  in  material  degrada¬ 
tion.  In  an  emergency  shutdown,  for  instance,  the  external 
load  is  disconnected  and  fuel  supply  shut  off.  As  long  as 
the  internal  thermodynamic  driving  potential  exists  to  sup¬ 
port  internal  reverse  current,  then  certain  portions  of  the  cell 
may  be  exposed  to  very  different  reactant  concentrations 
from  what  is  present  during  normal  operation.  For  exam¬ 
ple,  higher  oxygen  activity  may  be  produced  at  points  on  the 
fuel  electrode  having  positive  current  which  is  supporting 
the  negative  current  driven  elsewhere.  Clearly,  more  work  is 
needed  to  help  assess  the  true  impact  to  a  cell  from  these 
conditions. 

7.  Conclusions 

This  paper  examines  the  effect  of  load  transients  on  cells. 
Results  from  both  experiments  and  modeling  are  compared  to 
assess  the  ability  of  the  model  to  represent  fuel  cell  operation. 
Satisfactory  comparisons  are  found  between  the  data  and  the 
model.  The  model  is  then  applied  to  three  planar  geometries 
presently  being  considered  by  industry  (co-flow,  counter-flow 
and  cross-flow)  to  develop  improved  understanding  of  how 
SOFC  fuel  cells  respond  to  various  real-world  upset  condi¬ 
tions.  Of  particular  note  are  results  showing  load-decrease 
cases  having  reverse  current  over  significant  portions  of  the 
cell,  starting  from  the  moment  of  load  perturbation  up  to  the 
point  where  equilibrated  conditions  again  provide  positive 
current.  Such  reverse  current  conditions  arise  for  all  geome¬ 
try  types  studied,  although  the  cross-flow  case  is  more  severe 
for  the  same  magnitude  of  load  change.  An  estimate  of  the 
load  change  limit  such  that  negative  current  conditions  are 
avoided  is  80-85%.  Any  load  decrease  greater  than  this  will 
cause  current  reversal  over  at  least  portions  of  the  cell.  It  is 
suggested  that  conditions  for  current  reversal  be  examined 
further  given  the  potential  risk  for  electrode  degradation,  es¬ 
pecially  in  cases  of  emergency  shut-down  when  fuel  is  often 
removed  from  the  cell. 
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